A simple model of DNA denaturation and mutually avoiding walks statistics 
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Recently Garel, Monthus and Orland (Europhys. Lett. 55, 132 (2001)) considered a model of 
DNA denaturation in which excluded volume effects within each strand are neglected, while mutual 
avoidance is included. Using an approximate scheme they found a first order denaturation. We 
show that a first order transition for this model follows from exact results for the statistics of two 
mutually avoiding random walks, whose reunion exponent is c > 2, both in two and three dimensions. 
Analytical estimates of c due to the interactions with other denaturated loops, as well as numerical 
calculations, indicate that the transition is even sharper than in models where excluded volume 
effects are fully incorporated. The probability distribution of distances between homologous base 
pairs decays as a power law at the transition. 



PACS numbers: 87.14. Gg 05.70.Jk 05.70.Fh 87.15.Aa 



I. INTRODUCTION 

Simple models of DNA thermal denaturation have at- 
tracted a lot of attention for a long time jl]-f|] . Although 
it was pointed out rather early that excluded volume ef- 
fects may be crucial in determining the order of the de- 
naturation transition 0], it was only very recently that 
such effects were taken into account both at analytical 
and numerical level (5|-j^] . 

Even once accepted excluded volume effects as a key 
factor in determining the nature of the transition, one 
can still be interested in understanding how the relax- 
ation of self-avoidance constraints influences the overall 
behavior of the system. As this can be done in differ- 
ent ways, learning in detail about the effects of different 
geometrical exclusion constraints can help in conceiving 
better models or more efficient approximations. 

The main general aim of this paper is that of produc- 
ing a contribution to the debate on these general issues. 
In particular we analyze the model of DNA denatura- 
tion recently introduced and studied by Garel, Monthus 
and Orland (GMO) B. Among other things, we show 
that this model is a rather ideal context in which a re- 
cently proposed representation of DNA in terms of block 
copolymer networks Q can be applied. 

The Hamiltonian for the GMO model is M: 
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where fi(s), r 2 (s) describe the conformations of the two 
strands as functions of the curvilinear coordinate s and 
V(f f i2 (s)) is the attractive interaction between homolo- 
gous base pairs (ri 2 (s) = ri(s) — r 2 (s)). The last term 
on the r.h.s. of Eq. (Q) is due to excluded volume inter- 
actions, which act only between the two chains and not 



within a single chain. Thus the self-avoidance constraint 
is relaxed and each isolated strand follows a simple ran- 
dom walk (RW) statistics. 

The authors of Ref. || simplified further the Hamilto- 
nian (|l|) by approximating the excluded volume term by a 
long range interaction between the strands, which allows 
for an analytical solution of the problem. Within such 
approximation, the denaturation transition turns out to 
be of first order type j|, as in models where excluded 
volume effects are fully incorporated |s| ||. 

Our purpose here is to clarify further this point. We 
start from a description a la Poland - Sheraga (PS) |l]] of 
the GMO model and use a scries of exact results on re- 
union exponents of mutually avoiding walks, which show 
indeed that the denaturation transition is first order. 
More surprisingly, the transition turns out to be even 
sharper than that occurring in models where excluded 
volume effects are fully taken into account. This con- 
clusion is corroborated by approximate analytical treat- 
ments of excluded volume effects due to loop-loop inter- 
actions and by a scries of numerical results obtained for 
polymer lattice models in d = 2 and d = 3. 




FIG. 1. In the Poland - Sheraga model the full partition 
function of the DNA chain is factorized in terms of products 
of partition functions Sk of bound segments of length k and 
of partition functions Li of denaturated loops of length 21. 
Thick lines denote double stranded segments. 

In the Poland - Sheraga (PS) model the DNA chain 
is treated approximately as a sequence of non-interacting 
bound double segments and denaturated loops. Within 
this scheme the full partition function Zjf for a chain 
whose constituent strands have lengths N, is factorized 
in terms of elementary partition functions for loops and 
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segments (see Fig. Q). In the grand canonical ensemble, 
where a fugacity z per lattice step is assigned, the total 
partition function reads ||: 



= E 



- Zat — 



Vo(z) 



N 



l-L(z)S(z) 



(2) 



where -L(z) and «S(2) are the grand canonical partition 
functions for loops and double segments, respectively. 
The form of the numerator Vo(z) depends on boundary 
conditions and is not interesting for our purposes. The 
crucial quantity is the loop partition function, which in 
the most general form reads || : 



m = E 



(2iy 



(3) 



where /i is the (non-universal) effective connectivity con- 
stant, while c is a universal exponent. The order of 
the transition is related to the value of c J3|: At the 
transition temperature (which in the PS model corre- 
sponds to \j? z — > 1) loops are power-law distributed in 
size with probability P(l) ~ Z~ c , thus their average size 
(I) = ^2ilP(l) is finite for c > 2, implying a first order 
denaturation, while it diverges for 1 < c < 2 and the 
transition becomes continuous. 

For self- and mutually avoiding strands c was estimated 
analytically in an extended PS scheme ||, where ex- 
cluded volume effects between a self-avoiding loop and 
the rest of the chain are approximately taken into ac- 
count. This exponent was also obtained directly by 
Monte Carlo simulations |?],|| for polymer lattice mod- 
els. The agreement between the analytical estimates and 
simulation results for c is very good both in d = 2 and 
d = 3[@f. 

In the next Sections we will present analytical and nu- 
merical estimates of the exponent c for the GMO model. 
The former are based on the general theory of copolymer 
networks, which is briefly reviewed in the Appendix. 




FIG. 2. Example of loop in the GMO model. The two 
strands indicated as dashed and solid lines are allowed to 
overlap themselves, but are mutually avoiding. Equation (^) 
gives the associated reunion exponent. 



II. MUTUALLY AVOIDING WALKS 



with common origin and endpoint (see example of Fig. 
||). The probability that two random walks, with com- 
mon origin, meet for the first time after I steps decays 
asymptotically as P(l) ~ l~ c with 



c = dv — 1Vf}yy . 



(4) 



Here d is the dimensionality of the embedding space, 
v = 1/2 for RWs and rj[^ denotes the contribution of 
a vertex from which the two walks originate (see Fig. 
||) . Wc follow the notation of Ref. and indicate with 

( G) 

WfaL the exponent associated to a vertex from which /i 
solid and dashed lines depart, at the gaussian fixed 
point, which is the relevant one for our purposes. At this 
fixed point solid and dashed lines do not interact with 
themselves, but they avoid the lines of the other group. 

Figure ||, as well as Figs. |^-|| below, are examples of 
so called copolymer networks, i.e. networks whose con- 
stituents are polymers of different species. The relevance 
of these block copolymer networks for a description of de- 
naturating DNA has been already stressed in the context 
of models with full excluded volume effects |J. A gen- 
eral network of this type may be made, for instance, by 
an arbitrary mixture of random and self-avoiding walks 
(SAWs). Any pair of these walks may either avoid, or 
be allowed to cross each other. As for the homopolymer 
case, the general entropic exponent for a network with 
arbitrary topology follows from those of the constituents 
vertex exponents (see Refs. Jl0],[l3| for details). The gen- 
eral rule to calculate the entropic exponent of a network 
is simple and it is reviewed in the Appendix. Restrict- 
ing to the case in which all constituent polymers have 
the same v exponent, one associates to each independent 
loop a factor dv and to each vertex a factor —vq, with r\ 
the appropriate exponent which depends on the number 
of outgoing legs, and on their type. The generalization to 
an arbitrary mixture of polymer segments with distinct 
i/s, as for instance RWs and SAWs, is also possible. 

Coming back to Eq. (^), the vertex exponent for one 
solid and one dashed line is known exactly in d = 2 |L7) . 

(C) 

Its value is 77^' = —5/4, which leads to c = 2 + 1/4 > 2. 
In higher dimensions one has to resort to e expansion 
renormalization group results. In d — 3 resummation 
techniques yield f]^' « 0.57 which implies c w 2.07. 
Thus, both in d = 2 and 3 the reunion exponent is larger 
than the threshold value (c > 2), which immediately 
implies a first order denaturation for the GMO model, 
within the PS picture of non-interacting loops and bound 
segments. 

As a comparison we recall that for an isolated loop 
made by two self- and mutually avoiding strands the en- 
tropic exponent is [EH): 



dv. 



(5) 



A denaturated loop in the GMO model is formed by 
two mutually avoiding random walks of equal length /, 



For a SAW one has v = 3/4 in d = 2, which implies 
c = 3/2, while v ks 0.588 in d = 3, implying c w 1.76. 
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Thus, within the PS picture, denaturation is a continu- 
ous transition (c < 2) when full self- avoidance within a 
sine le loop is included 

It is also possible to go beyond the PS approximation 
in the GMO model, but it is already quite clear at this 
point that excluded volume effects between a single loop 
and the rest of the chain "localize" even more the loops 
and c can only increase. This can be shown explicitly for 
some simple geometries in which the loop is embedded. 
The calculations follow closely those of Ref. || for the 
model in which self - avoidance is fully included and are 
also briefly reviewed in the Appendix. We give here only 
the final results. 




FIG. 3. A small loop (between points A and B) embedded 
between two much larger loops is characterized by an expo- 
nent given by Eq. 

We consider the geometry shown in Fig. || in which 
a loop of length 21 is connected to two loops of length 
N — I each. In the limit I <C N the partition functions of 
the inner loop is still of the type of Eq. (||) and factor- 
izes with the partition function of the rest of the chain 
H . As explained in Appendix, the exponent c is however 
modified: 



c = dv — vrj^ . 



(6) 
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FIG. 4. In d = 2 there are two distinct vertices formed by 
two dashed and two solid lines (see text). 
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(7) 



Hence, for loops embedded in a geometry as in Fig. 
|3|, and formed by vertices of type (b) one finds c = 
dv — i/?74 MAW ^ = 3 + 5/8 = 3.625. Let us consider now 

vertices of type (a) for which we keep the notation ^2 
to indicate the associated exponent. Using the prescrip- 
tions of Ref. |17| one has 7722 = 35/12 and thus from 
Eq. (§) c = 2 + 11/24 « 2.46. Notice that stronger ex- 
cluded volume effects (as for vertices of type (b)) imply 
a higher value for c, with a quite large difference between 
the two estimated values. We point out that, obviously, 

(G) (MAW) . „ ,. . 

rjli — 772 m all dimensions. 



In d = 3 the above subtleties do not occur. The vertex 
exponent Tj^p « —1.81 was estimated in Ref. |l2j using 
e expansions and resummation techniques. Thus, from 
Eq. (0), one obtains c « 2.40. 

We consider now another geometry, where a loop is 
connected to two long double stranded segments, as il- 
lustrated in Fig. |H. Using the same method of Ref. 
we find (see Appendix): 



c = dv — 2vrj 12 



(8) 



where with rj 1 2 we indicate the exponent associated to 
a vertex where one double, and two single stranded seg- 
ments join. Notice that the double stranded segment is 
actually a self avoiding walk since, if it would overlap 
itself in some part, the mutual avoidance condition be- 
tween the two constituent strands would be violated. 




RW 



SAW 



FIG. 5. A loop embedded between two long double 
stranded segments. The latter follow the self-avoiding walks 
statistics. 



In d — 2 some care has to be taken in considering 
vertices, as the order of walks does matter. This is il- 
lustrated in Fig. ^ where the two possible distinct types 
of vertices with four outgoing lines, two solid and two 
dashed, are shown. For a vertex of type (a) more con- 
figurations are possible as neighboring lines of the same 
type are allowed to overlap. In (b), instead, the distinc- 
tion between solid and dashed lines becomes irrelevant 
as neighboring lines avoid each other. In the latter case 
the lines form a vertex of four mutually avoiding walks 
(MAWs) Q. A vertex with L outgoing MAWs, has an 
associated entropic exponent [^0 17 : 



Again in d = 2 one can use the techniques of Ref. fTi 



to find rj 1 



-2 - 7/16 which leads to c w 3.44. To 



our knowledge, there exists no field theoretical analysis 
for this type of vertices in higher dimensions, thus we 
are not able to make any predictions for c in d = 3, at 
present. 

Table | summarizes the exponents found in the GMO 
model for an isolated loop (c( lL >), for a loop interact- 
ing with other loops (c^ LL ^), as in the geometry of Fig. 
|3|, and for a loop interacting with double stranded seg- 
ments (c( LS )), as in Fig. ^. For the two dimensional case 
we reported both values of c^ LL - ) corresponding to the 
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loops formed by the two types of vertices of Fig. |4j. We 
notice that, as anticipated above, the interaction with 
other parts of the chain has the effect of increasing the 
loop exponent c^ LL \ c^ LS ) > c^ IL ^ as is the case for all 
other models studied so far The first order char- 

acter of the transition in the GMO model is therefore 
strengthened. 

It is interesting to compare the exponents obtained in 
this paper with the corresponding ones for the model in 
which strands are fully self-avoiding (FSA). The latter, 
taken from Ref. || , are given in Table The comparison 
between the two models indicates that the transition in 
the GMO model is generally sharper (higher c). This is 
always true with the exception of the exponent c( LL ' in 
d = 2 for vertices of the type of Fig. ||(a) for which we 
find c( LL > « 2.46, for the GMO model, while c< LL ) w 2.69 
in the FSA case. This is due to the fact that in the GMO 
model for the type of vertex considered excluded volume 
effects are not so pronounced as walks have considerable 
freedom to overlap. 

TABLE I. Summary of the exponents for the GMO model 
for an isolated loop (IL) and for a loop interacting with neigh- 
boring loops (LL) or segments (LS). The results for fully 
self-avoiding (FSA) strands (from Ref. ||) are also given. The 
last column shows the numerical estimates of c (those for the 
FSA are taken from Ref. ||). 





C (IL) 


C (LL) 


C (LS) 


c 


GMO (d=2) 


2.25 


2.46-3.62 


3.44 


2.95(5) 


FSA (d=2) 


1.50 


2.69 


2.41 


2.46(9) 


GMO (d=3) 


2.07 


2.40 




2.55(5) 


FSA (d=3) 


1.76 


2.22 


2.11 


2.18(6) 



III. NUMERICAL RESULTS 

We performed a series of numerical calculations for a 
lattice version of the GMO model both in d = 3 and 
d = 2 (cubic and square lattices). We considered two 
random walks of length N, described by the vectors f\{i) 
and r*2(i) (i — 0, 1, 2 . . . N) which identify the positions 
of the lattice monomers. The walks have common origin 
(ri(0) = ^2(0)), and no specific boundary condition is 
imposed on the other ends. While each strand can over- 
lap itself, mutual avoidance requires that r*i(i) ^ ^2 (J) 
for i =/= j. Overlaps at homologous points are however al- 
lowed (ri(i) — r*2 (i)); these correspond to a bound state 
of complementary base pairs, to which we assign an en- 
ergy e = — 1. At sufficiently low temperature (T) the 
walks are fully bound and form an object obeying self- 
avoiding walk statistics. 

We used the PERM (ll| algorithm and computed P(l) 
the probability distribution function (pdf) to find a loop 
of length 21 in the chain. While at low temperatures loops 
are exponentially distributed in size, at the critical point 



the pdf has an algebraic decay P(l) ~ Z — *\ from which we 
can extract the exponent c . Figure || shows a log- log 
plot P(l) vs. I in d = 3 at the estimated critical point 
(T c = 0.5181(3)), for chains of lengths up to N = 1280. 
A linear fit of the data yields c = 2.55(5), not far from 
the c (LL) given by Eq. (g) (see Table |). The result con- 
firms that the transition in the GMO model is sharper 
than when self-avoidance is fully included. Indeed, in the 
latter case the most accurate numerical estimate of the 
loop exponent is c = 2.18(6) g. 
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FIG. 6. Plot of log 10 P vs log 10 I at the estimated critical 
point in d — 3 and up to N = 1280. A linear fit of the data 
yields c = 2.55(5). 

We consider next P'(r), the probability that the dis- 
tance between homologous points in the strands \ri(k) — 
r<2,{k)\ equals r. This quantity was found to decay as a 
stretched exponential for r — > 00 || , from an analysis of 
the Hamiltonian of Eq. ([!]). 

This conclusion is at odd with the algebraic decay 
of the loop pdf found here for the GMO model, and 
also in other similar models with a first order transi- 
tion 0. For the numerical calculation of P'(r) we reg- 
istered all distances between homologous pairs along the 
chains \r\{k) — r2(fc)|, starting from k = 1 and up to 
the highest k for which the two strands are in contact 
(ri(fc) = T2{k)). Figure shows a plot of log 10 P'(r) vs. 
log 10 r at the transition temperature in d = 3. We cannot 
fit the data with a stretched exponential decay, but by 
increasing the system size the data approach a power law 
decay P'(r) - r~ k with k = 2.1(1). The fact that k > 1 
implies that the probability is normalizable. This is a 
confirmation that the transition is of first order type (see 
Ref. Notice however, that the algebraic nature of 

P'(r) implies infinite moments (r n ) starting from n = 2. 
We conclude that the exponential decay found in Ref. 
is an artifact of the approximation introduced. 

We performed a series of calculations also in the two 
dimensional case. At the estimated critical tempera- 
ture (T c = 0.5811(7)) we find, from the decay of the 
probability distributions P(l) and P'(r), the estimates 
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c = 2.95(5) and k = 3.0(1), confirming the first order 
nature of the transition. Again denaturation happens to 
be sharper than in the case where self-avoidance is fully 
incorporated, for which c « 2.46 g. In the PERM al- 
gorithm both types of vertices of Fig. |] are generated, 
thus c cannot be directly compared with any of the value 
obtained analytically in the previous section. 



o 




-10 1 1 1 1 1 1 1 1 1 

0.5 1 1.5 2 

lo Sio r 

FIG. 7. Plot of log 10 P' vs. log 10 r. The sets are the same 
as in Fig. ^. A linear fit yields k — 2.1(1). 

IV. CONCLUSION 

In this paper we studied the model for DNA denatura- 
tion recently introduced by Garel, Monthus and Orland 
H . Our analysis is complementary to that of Ref. || and 
it is based on the calculation of the exponent c describing 
the algebraic decay of the probability of denaturated loop 
lengths at the critical point. Using exact results from mu- 
tually avoiding walks statistics we find that c > 2 already 
for isolated loops, implying that the transition is of first 
order type. This conclusion is in agreement with the re- 
sults of Ref. p| , but it is not based on the approximate 
treatments pj| . 

For some specific simple geometries in which a loop is 
embedded in the chain we provide analytical estimates 
of c, following ideas from the theory of polymer networks 
HQ . These calculations show explicitly that c increases 
with respect to the isolated loop value, i.e. the transi- 
tion robustly maintains its first order character. More 
surprising is, at first sight, the fact that the denatura- 
tion in the GMO model is sharper than that occurring in 
models where excluded volume constraints are fully in- 
corporated (higher c) , as indicated by our analytical and 
numerical estimates of c both in d — 2 and 3. This be- 
havior can be explained by taking into account that the 
self- avoidance of the single strands somehow contrasts 
the excluded volume effects between each loop and the 
rest of the molecule. The result is that, when the single 
strand self-avoidance is relaxed, the loop sizes are further 
reduced, consistently with an increase of c. 



It is somehow remarkable how results of the theory of 
block copolymer networks, which are relevant to issues 
like that of establishing a link between field theory and 
multifractal measures (see e.g. Ref. |l^]), allow to draw 
reliable predictions on models of the denaturation tran- 
sition of DNA. 

There is another possibility of relaxing the self- 
avoidance which is somehow complementary to that of 
the GMO model. One can indeed neglect the mu- 
tual avoidance between the strands, while keeping self- 
avoidance within each stand. This possibility was stud- 
ied in d = 3 in Ref. and the denaturation transition 
turns out to be continuous (c < 2). Thus, relaxing the 
self- avoidance constraint in two different ways within the 
single loop, results in a softening (as in the latter ex- 
ample) or a sharpening (as in the GMO model) of the 
transition. 

Finally, we also considered the probability distribu- 
tion of the distances between homologous monomers of 
the two strands. We found that this quantity decays 
as a power law at the transition point, in disagreement 
with the stretched exponential behavior found in Ref. |j| , 
which is probably an artifact of the approximation intro- 
duced there. The algebraic decay of various quantities at 
denaturation, even in the case of first order transitions, 
is a common feature of other similar models. Finding 
the appropriate treatment of the Hamiltonian (|l|), which 
provides power-laws decay at the transition point is still 
an open question. 
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FIG. 8. Examples of polymer networks with (a) L — 2, 
ni = 5, 723 = 2 and 714 = 2, (b) C — 1, m = 2, 713 = 2 and (c) 
L = 3, 714 = 2. 
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APPENDIX 

We review here some known results for entropic expo- 
nents of networks of arbitrary topology. We focus mainly 
on networks whose constituent polymers are all of the 
same type, and just briefly mention copolymer networks. 

Let us consider a network made of self-avoiding walks, 
all avoiding each other, as the one shown in Fig. ||(a). Let 
N = J2i s i be the total length of the network, where each 
segment has length Sj. In the asymptotic limit N — » oo, 
Si —> oo, the total number of configurations for the net- 
work scales as fll3l : 



(9) 



where the factor fi is the effective connectivity constant 
for the walks (see e.g. Ref ]2(|) and / is a scaling func- 
tion. The universal exponent jg depends only on the 
number of independent loops, C, and the number of ver- 
tices with k outgoing segments nk as |Q 



7G 



1 - Cdv + ^2 v m n k 



(10) 



where rjk are exponents associated to a vertex with k 
outgoing legs. Such exponents are known exactly in 
d = 2 from conformal invariance [ fL3"[ , while they have 
been obtained from renormalization group and resum- 
mation techniques in d — 3 pit] . We notice that it is also 
customary to use the definition cr^. = vq^. 

This theory was recently applied to the denaturation 
problem |^|. In this case one considers network geome- 
tries as that of Fig. |](b) , in which a loop of total length 
21 is connected to two polymers of lengths N/2 each. 
Following the above notation one has C = 1, n\ = 2 and 
713 = 2 and thus the total number of configurations scale 



vertices with four outgoing legs, i.e. 77,4 = 2. In the limit 
I <C N one should find the partition function of two loops 
joined at a single vertex. Matching the scaling function 
in this limit yields 



(14) 



for a; — » 0. One finds thus that loops are distributed 
according to the probability P(l) ~ l~ c , with l|: 



c = dv — ur]4. 



(15) 



Finally, this theory can be generalized to the case of 
networks with more than one type of polymer, as for 
copolymer networks. Eq. (|l0|) can be generalized easily 
in the case that the constituting polymers have all the 
same metric exponents (for instance in the case of mutu- 
ally avoiding random walks) . The only difference is that 
now vertex exponents will not depend only on the num- 
ber of outgoing legs, but also on the type of polymers 
forming the vertex and, in two dimensions, also on the 
order. For instance, for a network made of two types of 
polymers, which are random walks and mutually avoid- 
ing with a geometry of Fig. |^ one arrives to a c exponent 
in the limit of an inner loop much shorter than the rest 
of the chain as given by Eq. ( 3|) , which is the copolymer 
network counterpart of Eq. (lq). Notice that while for 
homopolymers 772 = 0, in general a vertex with two out- 
going polymers of different type will have a non- vanishing 
exponent (as is the case of the 77^ introduced in Fig. ||). 

Finally if polymers constituting the inhomogeneous 
network will have different v exponent, as for networks 
formed by a mixture of random and self-avoiding walks 
than Eq. (^) is still valid, but in a form involving the 
radii of gyrations Ri of the constituent polymers, where 
Ri ~ s\, for a segment of length s». 



(11) 



with 7 = 1 — 6^ + 2^(771+773). Now in the limit I <C N one 
should recover the partition function for a self-avoiding 
walk of total length N which has an entropic exponent 
7 = 1 + 2^771. This requires that: 



x 



2l/T)3—dl/ 



(12) 



in the limit x — ► 0. Thus, loops of length I embedded 
within a long chain of length N I) are distributed 
according to the probability P(l) ~ l~ c , with [0]: 



c = dv — 2vrj3 



(13) 



As second example we consider a loop of length 21 em- 
bedded between two long loops of total length 2N, as 
depicted in Fig. ||(b). The total number of configura- 
tions is still given by Eq. ( pi] ) with 7=1 — 3dv + 21/774, 
as the network has now C = 3 independent loops and two 
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